function [AvE,Eps3,AvS,Sigma3] = ...
       avmatrix(coordinates,elements3,u,int_lambda,int_mu,...
               out_lambda,out_mu ,int_elements3,out_elements3);
%AVMATRIX  preprocessing for estimator and visualisation.
%   [AVE,EPS3,EPS4,AVS,SIGMA3,SIGMA4] = AVMATRIX(COORDINATES,ELEMENTS3,
%   ELEMENTS4,U,LAMBDA,MU) determines variables used as input for the
%   functions <aposteriori.m> and <show.m>.
%
%
%   See also FEM_LAME2D, AVMATRIX, and SHOW.

%    J. Alberty, C. Carstensen and S. A. Funken  07-03-00
%    File <avmatrix.m> in $(HOME)/acfk/fem_lame2d/cooks/ and
%                         $(HOME)/acfk/fem_lame2d/lshape_p1/ and
%                         $(HOME)/acfk/fem_lame2d/lshape_q1/ and
%                         $(HOME)/acfk/fem_lame2d/hole/

Eps3 = zeros(size(elements3,1),4);
Sigma3 = zeros(size(elements3,1),4);
AreaOmega = zeros(size(coordinates,1),1);
AvS = zeros(size(coordinates,1),4);
AvE = zeros(size(coordinates,1),4);

for j = 1:size(int_elements3,1)
  area3 = det([1,1,1;coordinates(int_elements3(j,:),:)'])/2;
  AreaOmega(int_elements3(j,:)) = AreaOmega(int_elements3(j,:)) +area3;
  PhiGrad = [1,1,1;coordinates(int_elements3(j,:),:)']\[zeros(1,2);eye(2)];
  U_Grad = u([1;1]*2*int_elements3(j,:)-[1;0]*[1,1,1])*PhiGrad;
  Eps3(j,:) = reshape((U_Grad+U_Grad')/2,1,4);
  Sigma3(j,:) = reshape(int_lambda*trace(U_Grad)*eye(2) ...
      +2*int_mu*(U_Grad+U_Grad')/2,1,4);
  AvE(int_elements3(j,:),:) = AvE(int_elements3(j,:),:) +area3*[1;1;1]*Eps3(j,:);
  AvS(int_elements3(j,:),:) = AvS(int_elements3(j,:),:) +area3*[1;1;1]*Sigma3(j,:);
end;

for j = 1:size(out_elements3,1)
  area3 = det([1,1,1;coordinates(out_elements3(j,:),:)'])/2;
  AreaOmega(out_elements3(j,:)) = AreaOmega(out_elements3(j,:)) +area3;
  PhiGrad = [1,1,1;coordinates(out_elements3(j,:),:)']\[zeros(1,2);eye(2)];
  U_Grad = u([1;1]*2*out_elements3(j,:)-[1;0]*[1,1,1])*PhiGrad;
  Eps3(j,:) = reshape((U_Grad+U_Grad')/2,1,4);
  Sigma3(j,:) = reshape(out_lambda*trace(U_Grad)*eye(2) ...
      +2*out_mu*(U_Grad+U_Grad')/2,1,4);
  AvE(out_elements3(j,:),:) = AvE(out_elements3(j,:),:) +area3*[1;1;1]*Eps3(j,:);
  AvS(out_elements3(j,:),:) = AvS(out_elements3(j,:),:) +area3*[1;1;1]*Sigma3(j,:);
end;

% for j = 1:size(elements3,1)
%   area3 = det([1,1,1;coordinates(elements3(j,:),:)'])/2;
%   AreaOmega(elements3(j,:)) = AreaOmega(elements3(j,:)) +area3;
%   PhiGrad = [1,1,1;coordinates(elements3(j,:),:)']\[zeros(1,2);eye(2)];
%   U_Grad = u([1;1]*2*elements3(j,:)-[1;0]*[1,1,1])*PhiGrad;
%   Eps3(j,:) = reshape((U_Grad+U_Grad')/2,1,4);
%   Sigma3(j,:) = reshape(lambda*trace(U_Grad)*eye(2) ...
%       +2*mu*(U_Grad+U_Grad')/2,1,4);
%   AvE(elements3(j,:),:) = AvE(elements3(j,:),:) +area3*[1;1;1]*Eps3(j,:);
%   AvS(elements3(j,:),:) = AvS(elements3(j,:),:) +area3*[1;1;1]*Sigma3(j,:);
% end;

AvE = AvE./(AreaOmega*[1,1,1,1]);
AvS = AvS./(AreaOmega*[1,1,1,1]);
